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Abstract. 

We make progress towards a 3D finite-element model for the magnetization of a 
high temperature superconductor (HTS): We suggest a method that takes into account 
demagnetisation effects and flux creep, while it neglects the effects associated with 
currents that are not perpendicular to the local magnetic induction. We consider 
samples that are subjected to a uniform magnetic field varying linearly with time. 
Their magnetization is calculated by means of a weak formulation in the magnetostatic 
approximation of the Maxwell equations (A-cj> formulation). An implicit method is used 
for the temporal resolution (Backward Euler scheme) and is solved in the open source 
solver GetDP. Picard iterations are used to deal with the power law conductivity of 
HTS. The finite element formulation is validated for an HTS tube with large pinning 
strength through the comparison with results obtained with other well-established 
methods. We show that carrying the calculations with a single time-step (as opposed 



Numerical simulation of the magnetization of high-temperature superconductors 2 
to many small time-steps) produce results with excellent accuracy in a drastically 
reduced simulation time. The numerical method is extended to the study of the 
trapped magnetization of cylinders that are drilled with different arrays of columnar 
holes arranged parallel to the cylinder axis. 
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1. Introduction 

Bulk high-temperature superconductors (HTS) become increasingly attractive for being 
used as efficient magnetic shields [I] or as powerful permanent magnets [2j [3]. Highly 
sensitive magnetic measurement systems, such as in a biomagnetic imaging device, 
need efficient magnetic shields for reducing the effects of the external magnetic 
disturbances [U [5j [6]. Powerful magnets are required in magnetic bearing systems, 
where they produce large levitation forces [7J |Sj, or in rotating machines, where they 
produce a large torque on the shaft [91 [TOj CD] • 

The performances of HTS trapped field magnets are limited by three main factors: 
(i) the critical current density in the sample, J c , that determines the maximum trapped 
magnetic field; (ii) the strength of the mechanical stresses, that arise from strong Lorentz 
forces and may result in cracks in the sample; and (iii) the heat exchange rate with the 
cryogenic fluid, that when too low may lead to significant temperature rises if the sample 
is subjected to a variable magnetic flux, as it is the case in rotating machines [T2] . 

In order to improve the performances of HTS magnets, attention has recently turned 
to bulk samples in which an array of parallel columnar holes is drilled along the c- 
axis [T3"j dU [T2] . The hole array enables one to obtain a better oxygen annealing [TB] 
- and therefore to raise J c — , to perform a more efficient cooling [17], or to reinforce 
mechanically the sample by injecting a resin [3j [15]. On the other hand, the holes also 
block the current stream lines — which have to flow around them — and, as a result, 
degrade the magnetic properties of the sample. In a previous work [18] , the Bean critical 
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state model was used for calculating the first magnetization of drilled cylinders of infinite 

extension, as a function of their hole pattern. It was shown that the penetration of the 

magnetic flux in a given hole generated a discontinuity in the flux distribution ahead of 

that hole. The trapped magnetic flux could then be increased by placing the holes on 

the discontinuity lines (i.e. lines where the current flow abruptly changes its direction 

due to the presence of a hole) of their direct neighbors, as this arrangement limited the 

perturbation by an individual hole of the overall flux distribution. 

This study was applied to samples of infinite extension and thus neglected 
demagnetisation effects. To pursue the study and compare with experimental results, 
work is needed to model the three-dimensional distribution of the magnetic flux while 
taking due account of the presence of the holes, the actual path followed by the current 
lines, and the resulting demagnetisation effects. 

Calculating a three-dimensional magnetic field distribution in HTS is notoriously 
difficult [121 [2U] • In the limit of infinite pinning strength, the magnetic flux distribution 
can be described with the concept of the critical state, introduced by Bean [2TJ. The 
critical state is characterized by a current density with a constant magnitude that flows 
perpendicular to the local flux density lines, since the magnetic force exerted on the 
vortices only depends on that component [20]. For a series of geometric configurations 
with a high level of symmetry, e.g. a cylinder subjected to an applied field with an 
axial symmetry, the current density is known to be everywhere perpendicular to the 
magnetic field, and the critical state can be simpled determined J2T], [22] ■ However, for 
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an arbitrary configuration where either the sample or the source of the field has no 
particular symmetry, the critical state model must be modified in order to properly 
describe the time evolution of the component of the current density that is parallel 
to the local magnetic field [20]. Moreover, for realistic pinning strengths, flux creep 
must also be taken into account, particularly at the temperature of liquid nitrogen 
(77 K) [23j [2U [25] . To our knowledge, no model includes yet a proper treatment of both 
the longitudinal component of the current density and the creep effects. 

The purpose of this work is to make progress toward a 3D model for the calculation 
of the trapped flux in drilled HTS magnets, by taking into account demagnetisation 
effects and flux creep, while neglecting the more delicate effects associated with 
longitudinal currents. In practice, we expect such a description to faithfully reproduce 
the actual flux distribution near the median plane of the sample, since the current lines 
are expected to lie in the plane and thus to be perpendicular to the local flux lines. 
Moreover, we believe that the neglect of those effects associated with the longitudinal 
component of the current density leads nonetheless to a first-order approximation from 
which qualitative conclusions can be drawn regarding the influence of the hole lattice 
on the magnetic properties of the drilled samples. 

In the rest of this paper, HTS are modeled by a power law conductivity cr(E) [25] 
assuming the form 



where E c is the critical electric field and J c is the critical current density. The critical 
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exponent, n, is related to the pinning strength in the material and is assumed to be 

independent of the magnetic field. When n = 1, we recover the constitutive law of 

an ohmic material. In the opposite limit of infinite pinning strength, n — > oo, the 

power law model is asymptotically equivalent to the Bean model [2SI [27]. Several 

numerical methods are available to solve for the magnetic field penetration in HTS 

with a power law conductivity: finite-difference approximation in cylinders [28] , Green's 

function approach in cylinders [IH] or in tubes pp, and finite element method (FEM) 

with so-called A-(j) formulations [2H1 ED], formulations [3TJ 132] , or unconstrained 

//-formulations [33l El] . In each of these methods, the choice of the time-step is crucial 

since it governs the convergence rate and the total calculation time, which can become 

excessively long on a 3D mesh when n is large [301 ESS 135] . 

To our knowledge, in the FEM suggested so far, the computation time-step was 

chosen much smaller than the timescale characterising the simulated external excitation. 

Such a choice can however be largely improved in the (present) case of an excitation 

varying linearly with time. Our argumentation is two-fold. First, from the point of view 

of the physics involved, one knows that the vortex motion and flux creep are strongly 

reduced as the pinning strength increases. Thus, for large n, the motion of vortices can 

only be induced by applying an external flux variation, so that the time behaviour of the 

magnetic response is expected to be mainly dictated by the excitation sweep rate, not by 

creep effects. The second part of our argumentation stems from the numerics involved. 

We solve a time-differential equation of the form du/dt = g(u) with the backward Euler 
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scheme |36j. The temporal derivative at time t is approximated at first-order, yielding 

the implicit equation 

At =9{Ut) - (2) 
Such a scheme has been shown to yield a truncature error proportional to the second time 
derivative, e t « d 2 u/dt 2 At + 0(At 2 ) [36]. A gain, in the limit of large pinning strength 
and with an external field applied as a ramp, we expect the second time derivative of 
the magnetic response to be small, as its timescale is dictated by that of the excitation, 
which varies here linearly with time. These arguments suggest that larger At can be 
used, and in the extreme case, a single time-step might be used. 

This paper addresses the questions of the accuracy and the convergence of a 
single time-step method that is suitable for a 3D model of HTS. For this purpose, 
we use a finite-element formulation implemented in the open source numerical solver 
GetDP [371 [38] . The rest of this manuscript is organised as follows: in Section El we 
describe and motivate the choice of an A-<f> formulation. In Section [31 we describe the 
implementation of this formulation into GetDP, and validate it in Section [H where 
comparisons are made with the Bean model in the case of an HTS tube with an infinite 
height (2D geometry), and with the Green's function method [19] in the cases of a tube 
of finite height (3D geometry). In particular, we analyse the validity of the single time- 
step method as a function of the value of the critical exponent n and the ramping rate. 
In section [SJ we apply the FEM for calculating the trapped magnetic flux density in 
drilled HTS cylinders with a finite height, for four different periodical arrangements of 
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the columnar holes. We then conclude in Section 



2. Finite element A — (ft formulation 

The description of the magnetic field penetration in HTS is based on magneto-quasistatic 
approximation of the Maxwell equations [39]. The HTS conductivity is given by 
Equation (JTJ) and the lower critical field, H&, is neglected against the applied field, 
so that the material follows the constitutive law, B = /ioH. We introduce the vector 
potential A and the scalar potential <f>, through 

B = B sclf + 5 a (t)e z = Vx A + Vx A a , (3) 

E=-^-^-V* (4) 

dt dt r ' v ; 

where the magnetic flux density is split into two contributions: the uniform applied 
magnetic flux density, B a (t)e z , which points along the z-axis and varies linearly 
with time as B a (t) = B a t, and the reaction magnetic flux density, B se if, which is 
produced by the eddy currents induced in the HTS. In cylindrical coordinates, the 
vector potential corresponding to the uniform applied magnetic flux density is given 
by A a = —r/2B a (t)ee and we have dA a /dt = —r/2B a ee. The introduction of the 
potentials A and <fi into the magneto-quasistatic Maxwell equations leads to two coupled 
equations: 

V x V x A = /i c(A, 0) (—A - A a - V0) , (5) 



V • {a(A, 0) (-A - A a - W) } = 0, (6) 
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where the electrical conductivity o is calculated from the power law (pQ), as a = 

J e /EW n '*\E\V-- n V n . These equations are sufficient to describe the electromagnetic 

behavior of HTS in the A-<f> formulation [291 SO]- The choice of this particular 

formulation is motivated by the fact that it produces a strong knowledge of the magnetic 

flux density, which is the quantity that is directly available in experiment. The Dirichlet 

boundary conditions on A and </> are imposed on the outer surface of a circular shell 

(in 2D geometry) or a spherical shell (in 3D geometry), whose external surface in both 

cases is sent to infinity by a Jacobian transformation [41 J. At infinity, we set A = and 

= 0. 

Equations ([MS]) are solved by the Galerkin residual minimization method, which 
yields 

(V x A, V x A { )- < B self x n, A, > -/i (a(A + A a + V0), A;) = 0, (7) 

(a A, V0j-) + (aA a , Vfa) + (aV0, V^)- < n.E, fa >= 0, (8) 

where and <f>j are basis functions that are known a priori, the notation (u, v) 
corresponds to the volume integral f n uv dV over the volume Q, and < u, v > stands for 
the surface integral f dn uv dC. Surface terms are used for imposing Neumann boundary 
conditions when appropriate. 

The vector potential is approximated as a series of basis functions, A = J2i ^iAj, 
where the Aj's are first-order edge functions that ensure the continuity of the normal 
component of B from mesh to mesh. The vector potential is defined in a gauge for 
the edge functions that avoids the computation of the curl of A for the post-processing 
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of the magnetic flux density [42] . This gauge condition reads A.w = 0, where w is a 

set of edges that connects all the nodes of the mesh through an open path, in such a 

way that a given pair of nodes can only be connected by a unique continuous path (see 

Figured]). Similarly to A, the scalar potential is expanded as (f> = J2j bj<f>j, where fa 

are first-order nodal functions. 

3. Computation of the finite element model 

The weak formulation ([318]) is implemented into the open-source solver for discrete 
problems, GetDP [37, 38J. GetDP presents two major advantages over commercial 
finite-element softwares: it is available free of charge and it offers a large choice of 
numerical methods to be implemented with a full control of the inherent parameters. 

As stated in the Introduction, we use a step-by-step temporal resolution with a 
backward Euler scheme, which has a good stability and a high convergence rate even 
with very large time steps (36] . The convergence and the stability of this method has 
already been demonstrated in the context of HTS in the case of a E'-formulation |43j. 
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Figure 1. Example of a set of meshing edges w used for the definition of the vector 
potential A. 
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Figure 2. Operations to be executed for one time iteration. 

In our formulation, the implicit resolution required at each step generates a system of 
equations which are non-linear, because of the conductivity law of Equation ([ID. This 
non-linearity is treated with a Picard iterative loop JH], which consists in updating at 
each time step the value of the non-linear term with the solution found at the previous 
iteration. The loop is run until the relative difference between two consecutive solutions, 
e n , is smaller than a predefined criterion, taken empirically here as e n < 2.10 -3 . Using 
a Picard iteration scheme with a power law conductivity prevents one from having to 
deal with the infinite derivatives that appear in the more traditional Newton-Raphson 
scheme |45J. Figure [2] schematically represents the sequence of operations to be executed 
during a given time step. 

In the following, we will compare two different choices for the time integration of 
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the field from instant "zero" to a predetermined instant, t\: in the first choice, the 

integration is carried out in a succession of small time steps of duration At <C t±; in the 

second choice, the equations are iterated in a single time-step, with At — t\. These two 

choices will be compared in a number of situations. 

4. Simulation of the magnetization of a HTS tube 

We first apply the FEM to the calculation of the magnetization of an HTS tube subjected 
to an axial field, in both limits of infinite and finite height. These are geometries for 
which the current density is everywhere perpendicular to the local magnetic field and 
solutions are known from other methods. The goal of this section is to compare the 
FEM to these other methods to validate our approach. 

The high level of symmetry in each geometry allows us in principle to reduce the 
mesh dimension. However, we deliberately choose not to exploit symmetry to construct 
the mesh, so as to use the weak formulation ((7j)-(jHJ) without simplification since that 
formulation will be used in geometries having no such symmetries. Thus, for the case 
of a tube of infinite extension, we use a 2D mesh of the cross section, while for the case 
of a tube with a finite height, we use a 3D mesh. The FEM results are compared to the 
predictions of the Bean model in the case of infinitely long tubes and to the results of 
the Green's function of Brandt [19] in the case of tubes with a finite height. 
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4-1. Tube of infinite extension (2D geometry) 

Consider first a superconducting tube of infinite height subjected to a uniform magnetic 
field applied parallel to its axis, as a ramp. The tube has an external radius a = 10 mm 
and an internal radius 6 = 5 mm. The pinning forces in the superconductor are assumed 
to be infinite. Under this hypothesis, the Bean model [21] applies and predicts that the 
magnetic flux density decreases linearly inside the wall and is constant inside the hole. 

As explained in Section [lj the power law conductivity jT]) is asymptotically 
equivalent to the Bean model when n — > 00. From a practical point of view, it has been 
shown in Ref. [IB] that the use of the power law with a critical exponent of n — 100 and 
with a sweep rate of B a = 10 mT/s yields an accurate approximation of the Bean model. 
In this section, we choose that parameter values in order to compare the results of the 
finite-element model to analytical expressions of the Bean critical-state. The critical 
current density J c is assumed to be independent on the magnetic flux density! and has 
a value of J c = 2 10 7 A/m 2 . The critical electric field E c is taken to be E c = 10 -4 V/m. 
The theoretical penetration field of the tube H p , is given by H p = J c (b — a) = 10 5 A/m, 
which corresponds to a flux density, [IqH p = 125.6 mT. 

In Figure [3]- (a), the magnetic field profile is plotted along the diameter of the tube 

for an external induction B a = 10, 50, 100, 150, and 200 mT. FEM simulations are run 

with different choices of time-steps: dashed lines show the results of simulations with 

multiple small time-steps At = 1 s, stopping at either ti = 1, 5, 10, 15 and 20 s; solid 
f Note that a model with field-dependent J c can easily be implemented in GetDP 
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Figure 3. (a) Profile of the magnetic flux density along the tube diameter, as 
calculated with the FEM with a single time-step (solid lines) and with multiple time- 
steps At = 1 s (dashed lines). The critical exponent n is chosen large (here n — 100) so 
as to approach the critical state. The magnetic flux density is applied with a constant 
sweep rate of 10 mT/s. The profiles are shown for B a = 10, 50, 100, 150, and 200 mT. 
(b) Average deviation AB from the Bean model, as a function of B a , for the single 
time-step method (circles) and for the multiple time-step method (squares). 

lines show results from single time-step simulations, where At = t\ is fixed to either 1, 
5, 10, 15, or 20 s. It can be observed that in each case the profile of the magnetic field 
in the superconductor is linear. It closely follows the result of the Bean model, 

^Bean = B a - ^ J c {a - t) (b < T < a) , 

(9) 

^Bean = B a - fl J c (a - b) (0 < T < b) . 

To further quantify the results, we define the average deviation from the Bean model as 

1 ra 



AB 



2a 



\B 



FEM 



BBca,n\ dr, 



(10) 



where -Bfem stands for the FEM results. Figure [3]- (b) shows the average deviation in 
the FEM method using a single time-step (filled circles) and in that using multiple time- 
steps (open squares). Both methods produce almost the same deviation as long as the 
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magnetic field has not fully penetrated the wall of the tube, or for B a < 120 mT. For 

B a > 120 mT, the error obtained with the multiple time-step approach first increases 

abruptly and then has a value around 3 mT. The single time-step method, on the other 

hand, leads to an error which peaks at about 2.5 mT at B a = 130 mT and then decreases 

at larger fields to be less than 2 mT at B a = 200 mT. 

Overall, the single time-step method gives an accurate solution with a relative 

error that stay below 5% of the Bean prediction. Note that this result is obtained in 

a relatively short calculation time with respect to a multiple time-step method. For 

example, the 20 simulations of Figure [3]-(b) take less than half a day on a dual-core 

2.8 GHz processor with 2 Gb of memory, whereas the multiple time-step approach with 

At = 1 s takes almost 3 days on the same computer. 

4-2. Tube of finite extension (3D geometry) 

We now turn to the case of a superconducting tube of finite height subjected to a uniform 
axial field. The tube has an external radius of 10 mm, an internal radius of 5 mm (see 
Figure HJ), and a height of 8 mm. The external field is applied with a constant rate 
B a = 10 mT/s and raises up to B a = 200 mT. Here again, a large pinning strength 
with n = 100 is assumed. The critical current density, J c and the critical electric field, 
E c have the same values as for the tube with an infinite height. 

The FEM approach is carried on a 3D mesh with a single time-step method 
(At = 20 s). Only half of the tube is actually meshed, and vanishing conditions on 
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Figure 4. Sketch of the HTS tube. The outer radius is 10 mm, the inner radius is 
5 mm, and the height is 8 mm. The arrows indicate different scan directions for plotting 
the magnetic flux profile. The external field is applied with a constant sweep rate of 
B a = 10 mT/s, and raises up to B a — 200 mT. The flux penetration problem is either 
solved with the FEM single time-step method (gray solid lines and At = 20 s), or with 
Brandt's method with multiple time-steps (black dashed lines and At — 5 10~ 4 s), 
both assuming n = 100. 

the tangential component of B are imposed in the median plane. The FEM results 
are compared with those of the Green's function method of Brandt [TJ [19], which in 
this geometry is based on a 2D-kernel. The time-step is fixed at 5 10~ 4 s to ensure 
convergence. 

The z-component of the magnetic flux density is probed along three different 
directions [see Figure HJ (a), the tube axis, (b) a diameter at the top surface, and 
(c), a diameter on the median plane]. Solid lines show the FEM results and dashed lines 
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show those of Brandt's method. It can be observed that solid lines exhibit step-like 

features, whereas dashed lines are smoother. This difference can be traced back to the 

low meshing density adopted in the FEM method. Even though only half of the tube 

is meshed in 3D, the maximum number of available nodes with 6 Gb RAM (200 000) 

is still too small to obtain a smooth curve after linear interpolation, unlike the Brandt 

method which uses a specific interpolation on the 2D meshing. 

Despite this observed difference, one can see that the results of the two methods 
are in good agreement. In particular, on the linescan along direction (c) (Figure ffl(c)), 
we observe a magnetic flux density on the outer wall of the tube that is slightly larger 
than B a = 200 mT. This is caused by the demagnetizing field, which was absent in the 
results of the previous subsection. Inside the superconducting wall, the magnetic flux 
density decreases linearly; in the central part of the tube, it remains at a low level, but 
exhibits variations due again to the demagnetizing field. 

Figure O shows the z-component of the magnetic flux density calculated at the center 
of the cylinder, S cen t er , as a function of the external field B a . The dashed lines show 
the result of the Green's function approach. Circles show the FEM results in a single 
time-step approach, with different choices of the time-step ranging between At = 1 s 
and 20 s. Here again, the agreement between the methods is excellent, demonstrating 
the relevance for adopting a single time-step iteration in a FEM approach. 
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Figure 5. Magnetic flux density calculated at the center of a HTS tube with the 
FEM single time-step method (circles) and with the Green's function method (dashed 
lines and At = 5 10~ 4 s). The external field is applied with a ramp of constant rate 
B a = 10 mT/s and increases up to 200 mT. FEM single time-step method are shown 
for different choices of the time-step: At — 1, . . . , 20 s. 
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Figure 6. Difference between the magnetic flux density at the center of the HTS tube, 
as calculated by the FEM single time-step method and the Green's function method 
(with multiple time steps At = 5 10 -4 s). The magnetic field is applied as a ramp with 
a sweep rate of 1 mT/s (squares), 10 mT/s (circles) and 100 mT/s (triangles). The 
applied magnetic flux density is ramped up to B a = 200 mT. The critical exponent n 
varies from 5 to 100. 
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4-3. Domain of validity of the single time-step method 

We have seen in the two previous subsections that the FEM method with a single time- 
step produces accurate results in the strong pinning limit (n = 100). The purpose of 
this subsection is to analyze the accuracy of the method at lower pinning strength and 
establish its sensitivity to the sweep rate of the external field. Mastering these two 
factors is essential to make comparisons with experiments. 

We estimate the error of the FEM single time-step method on the basis of the 
magnetic field produced at the center of the tube with a finite height (the tube considered 
in the previous subsection). The external field is ramped with a fixed rate B a up to 
B a = 200 mT. The error is then evaluated as the absolute difference between the results 
of the FEM and the Green's function methods. Figure [6] shows the error (in %) as a 
function of the critical exponent n, varying from 5 to 100, and the sweep rate B a , taken 
as 1 niT/s (circles), 10 mT/s (squares), or 100 mT/s (triangles). 

In the very strong pinning limit (n = 100), the error remain small (below 2 %) and 
is fairly independent of the sweep rate. This limit corresponds to the critical state which 
is uniquely determined by the external conditions and is independent of the magnetic 
history of the sample. Provided convergence is guaranteed, the FEM approach should 
thus produce the critical state solution. The opposite limit of low pinning strength shows 
a much larger sensistity to the sweep rate and a larger spread in the error. Here, these 
results should be considered as qualitative only, as the Green's function method itself 
has an error that grows in this limit, so that our estimate of the FEM error becomes 
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questionable in this regime. For intermediate values of n, the error remains low and 

weakly sensitive to the sweep rate: e.g., for the experimentally relevant value for melt- 
textured YBCO at 77 K, n = 20 [17], the error is below 3 %. This demonstrates that the 
single time-step method is useful for simulating the magnetization of HTS with finite 
pinning strengths. 

5. Magnetization of drilled cylinders 

An extension of the single time-step method is presented in this last section where we 
compare the magnetization of cylinders containing 4 different arrays of holes. In a 
previous work [18], the Bean critical state has been used to compare the magnetization 
of cylinders of infinite height with four different patterns of holes: the squared and the 
centered rectangular lattices having a translational symmetry, and the polar squared 
and polar triangular lattices with a rotational symmetry. It was found that the largest 
trapped magnetic flux is obtained with the polar triangular lattice. We now consider 
FEM calculations in order to take into account demagnetisation and creep effects. 

To this end, we consider cylinders (radius of 10 mm and height of 8 mm) that are 
drilled by the four lattices considered in Ref. [IB]. The lattice parameters are chosen 
in a such way that the total diameter of the holes is constant (507r mm), so as to fix 
the total surface of heat exchange. The squared and the centered rectangular lattices 
contain each 25 holes with a radius of 1 mm. The polar lattices contain two layers of 10 
holes with a radius of 1 mm and a central layer with 10 holes with a radius of 0.5 mm. 
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The four samples are represented in Figure [3- (a). 

In order to calculate the trapped magnetic flux, HTS samples are first magnetized 
by an external field varying linearly with time. A magnetic flux is then trapped in the 
sample when the external magnetic field returns to mT. This magnetization process is 
calculated here in two time-steps: one for increasing the applied magnetic flux density 
to 600 mT with a constant sweep rate of 10 mT/s and a second one for decreasing it to 
mT with the same sweep rate. 

Figure EJ-(b) shows the trapped magnetic flux density profile along the cylinder 
diameter [see black arrow in Figure 0-(a)]. The dashed curve corresponds to the 
trapped flux profile in a cylinder having the same geometry and material parameters, 
but containing no holes. The flux profiles exhibit steps, resulting from the low number 
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Figure 7. (a)- The four lattices in the cylinder (radius 10 mm, height 8 mm). (b),(c) 
- Trapped magnetic flux profile along the cylinder diameter, as calculated with the 
FEM two time-steps method, with a sweep rate of 10 mT/s and n = 25. The profile is 
represented in the top cross section of the cylinder. The dashed line corresponds to the 
flux profile in a cylinder without holes, also calculated with the FEM two time-step 
method. 
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fimax (mT) 


3D - top 


3D - CENTER 


Infinite height 


Polar triangular lattice 


70.05 


112.7 


137.9 


Polar squared lattice 


63.6 


97.7 


120.9 


Squared lattice 


56.7 


87.8 


110.8 


Centered rectangular lattice 


50.6 


76.4 


101 



Table 1. Maximum trapped magnetic flux density in cylinders of finite height, as 
calculated in the top cross section (3D - top) and in the median plane (3D - center), 
and in cylinders of infinite height obtained in Ref. |18j . 

of meshing elements used in 3D simulations as was already observed in Section 14.21 
It can be observed that the maximum trapped magnetic flux density is smaller in the 
drilled samples than in the bulk one. 

Table [1] lists the values of the maximum trapped magnetic flux density in the 
top cross section and in the median plane, as well as the results obtained for infinite 
cylinders |18j . In all cases, the maximum trapped magnetic flux density is obtained with 
a polar triangular lattice, with a value higher by ~ 40% with respect to that obtained 
in a centered rectangular lattice. This result is independent of the cross-section where it 
is calculated and agrees with the theoretical predictions based on the Bean model [18J. 
The demagnetisation effects only affect the values of the maximum trapped flux density 
that are smaller in the finite height samples than in the cylinders of infinite height with 
the same hole lattice. 
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6. Conclusions 

Using the open source solver GetDP, we have implemented a 3D finite element A-(f) 
formulation for the calculation of the magnetization of bulk HTS subjected to a ramp 
of magnetic field. The numerical method is based on a single time-step iteration that 
reduces drastically the total calculation time. By comparing it to the Bean model [2TJ 
in infinite tubes and to the Green's function method [19] in tubes of finite height, we 
have shown that the FEM approach accurately describes the magnetic properties of 
superconductors with strong pinning. Although it neglects the effects associated with 
currents that are parallel to the magnetic field, that study makes progress toward a 3D 
model of HTS magnets that takes into account demagnetisation effects and flux creep. 

As an extension of the FEM single time-step method, we have calculated the 
trapped magnetic flux in drilled cylinders of finite height. The numerical method uses 
only two time-steps: the first one during the ramping up of the applied field to i? max 
and the second one for the return of the external field to zero. Using this method, we 
have been able to extend a previous analysis for tubes of infinite extension to a full 3D 
geometry. These results confirm that the trapped magnetic flux is maximized by drilling 
the holes according to a polar triangular lattice. 
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